Spontaneous Crystallization of Skyrmions and Fractional Vortices in the Fast-rotating 
and Rapidly-quenched Spin-1 Bose-Einstein Condensates 
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We investigate the spontaneous generation of crystallized topological defects via the combining 
effects of fast rotation and rapid thermal quench on the spin-1 Bose-Einstein condensates. By 
solving the stochastic projected Gross-Pitaevskii equation, we show that, when the system reaches 
equilibrium, a hexagonal lattice of skyrmions, and a square lattice of half-quantized vortices can be 
formed in a ferromagnetic and antiferromagnetic spinor BEG, respetively, which can be imaged by 
using the polarization-dependent phase-contrast method. 
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Topological defects are a manifestation of sponta- 
neously broken symmetries [1 . Formation and obser- 
vation of topological defects are one of the most fun- 
damental and fascinating topics in various aspects of 
physics, ranging from condensed matter physics to cos- 
mology. However, owing to the limitation of energy scales 
in the earth-bound physics experiments, topological de- 
fects are mostly created and observed in the condensed 
matter systems. For example, magnetic domains walls of 
magnetized material and string defects in ^He superfluid 
phase transitions have been extensively studied [2 . 

Recently, owing to the realization of spinor Bose- 
Einstein condensate (BEG) of alkali atoms in optical 
trap [21 [4], the creation of topological defects in ultra- 
cold atomic systems has become possible. A spinor BEG 
is fully characterized by the spin degrees of freedom and 
behaves as a vector in the spin space. Theoretical stud- 
ies for the spinor BEG were pioneered by Ho [5 , and 
independently by Ohmi and Machida [6 . More recently, 
spinor BEGs with F > 1 [7H3 or with long range dipolar 
interaction [TOl [11] have also been theoretically investi- 
gated. In general, the physical behavior of the spinor 
BEG depends crucially on its magnetic properties, such 
that the interplay between the superfluidity and mag- 
netism of the condensate makes the spinor BEG a can- 
didate to exhibit a variety of nontrivial ordered states, 
such as the skyrmions [12], Mermin-Ho vortices [13], 
Anderson-Thouless vortices [M] and so forth. 

So far, all theoretical studies regarding the formation 
of topological defects in spinor BEG appeal to manipu- 
late the external and internal degrees of freedom of the 
condensed atoms at zero temperature. On the other 
hand, according to the Kibble-Zurek scenario, topolog- 
ical defects can also be created through phase transitions 
at finite temperatures, which are fundamentally caused 
by spontaneous symmetry breaking and thermal fluctua- 
tions near the critical point. In this paper, we show that 
it is possible to create crystalline orders of skyrmions 



and fractional vortices simply by thermally quenching a 
rotating spin-1 BEG. This enables us to probe into the 
very fundamental aspects of topological defects without 
any engineering of dynamical processes, since evaporative 
cooling is prerequisite in creating BEGs and the methods 
of rotating condensates have been well developed in a va- 
riety of ultracold atoms experiments. In the framework of 
mean field theory, the dynamics of a BEG at nonzero tem- 
peratures can be described by the stochastic projected 
Gross-Pitaevskii equation (SPGPE) [15], which relies on 
the assumption that the system can be treated as a con- 
densate band in contact with a thermal reservoir compris- 
ing of all non-condensed particles. In such a scheme, the 
condensate band is described by the truncated Wigner 
method [16 including the projected c-field method, while 
the non-condensate band is described by the quantum ki- 
netic theory [T7l[T8]. In the following, we shall solve the 
SPGPE numerically for a rotating trapped spin-1 BEG. 
We show that when the system is quenched down to a 
very low temperature, a lattice of skyrmions and half- 
quantized vortices (HQVs) can be created in the spinor 
BEG of ^^Rb and ^^Na, respectively. 

The spin-1 BEG is characterized by a vectorial or- 
der parameter, ^ = ( ^i, ^o, ^-i ) (the superscript 
T stands for the transpose), where {j = ±1,0) denotes 
the macroscopic wave function of the atoms condensed in 
the spin states, \F = l,mi? = ±1,0), respectively. The 
dynamics of in a confining potential is described by 
the following coupled nonlinear Schrodinger equations 

ihdt^i = HfP^i (1) 

a=x,y,z n,k,l=0,±l 

where ?4 = —h^^'^/2m-\-V{T)-\- Qn |^|^ denotes the spin- 
independent part of the Hamiltonian, and Fq, are the ma- 
trices representing the Gartesian components of the spin 
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angular momentum with quantization axis fixed in the 
2;-axis. The couphng constants Qn and Qs characterizing 
the density-density and spin-exchange interactions, re- 
spectively, are related to the 5-wave scattering lengths ao 
and a2 in the total spin channels Ftotai = and Ffotai = 2 
by gn = 47rfi^ (ao + 2a2) /3m, gs = 4:7rh^ (a2 — ao) /3m 
[5j [6j. The ground state of the spinor BEC depends 
crucially on the sign of gs. For the ferromagnetic cou- 
pling, gs < 0, the condensate is in the "axial" state. 
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= 1 [5 . On the other hand, for the antiferro- 
magnetic coupling, > 0, the condensate is in the "po- 
lar" state, ^F^ =0 [5 . As we shall focus on the vor- 
tices formed by the spin textures, it is more convenient 
to introduce the basis set = ^^V^^)^ such that 

^±1 = + /V^ and ^0 = "^z- As a result, 

Fa\<^) = 0, and the spin texture, which is parallel to the 
local magnetic moment, is defined by S (r) = ip~^^^ 
where ^=(^3^,^2/5^2) [19- Consequently, we have 

and Sz oc j^*! |^ — |^. For later convenience, we define 
the unit vector s (r) = S (r) / |S (r)|. 
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FIG. 2: (Color online) (a) The equilibrium spin textures for 
the spinor BEC of ^^Rb; (b) Distribution of in the spin 
textures; (c) The orientations of the unit vector s (r) for a 
skyrmion. The color of each arrow indicates the magnitude 
of Sz\ (d) The density profile of the topological charges. 
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FIG. 1: (Color online) Snapshots of (a) 
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(c) |^i|^; (d) 1^1 for the quenched rotating spinor BEC of 
^^Rb. When the system reaches equilibrium at T = 10 nK 
with Q — 0.95, /i = 8 (rightmost column), crystalline order of 
vortex-trimers is established in each ^ j . The particle numbers 
in the spinor BEC are iV±i 4.84 x 10^, A^o ~ 4.79 x 10^. 



To study the non-equilibrium dynamics of a quenched 
rotating spinor BEC, we shall generalize the formulae in 
Ref.[15^ to the following set of coupled SPGPEs 



(2) 



a ^ 



7j 



dt 



} 



where T and /i denote the final temperature and chemi- 
cal potential, jj the growth rate for the j-th component. 



and dWj/dt is the complex- valued white noise associ- 
ated with the condensate growth. The projection op- 
erator V restricts the dynamics of the spinor BEC in 
the lower energy region below the cutoff energy E^. In 
the rotating frame, ?^ is replaced by ?^ — 1^1/^, where 
Lz = —ih {xdy — ydx) is the z-component of the orbital 
angular momentum, and Vt is the angular frequency of 
rotation. As we shall focus on the fast rotating BECs, in 
which the atomic cloud appears highly pancake-shaped, 
it is reasonable to treat the system as two-dimensional. 
We therefore assume V(y) = muj'^ (x^ -\- y'^ + A^2:^) /2 
with A ^ 1. The effective 2D interaction strength can 
be obtained by integrating the wave functions with re- 
spect to z to eliminate the axial degree of freedom. The 
numerical procedures for integrating the set of coupled 
SPGPEs, are described as follows. First, the initial state 
of each is sampled by using the grand-canonical en- 
semble for free ideal Bose gas at a temperature To be- 
low the critical temperature and of chemical potentials 
/jjj^Q. The spatial dependence of the initial state can be 
specified as a linear combination of some basis functions. 
Here, we adopt the basis consisted of plane waves with 
quantized momentum k = 27r {nx^riy) /L (n^^, Uy are in- 
tegers and L is the size of the computation domain), i.e., 
{t = 0) = Xlk^ %;k'0k (r), where ipi^ (r) are the plane- 
wave basis functions. The condensate band lies below the 

an energy cutoff Er > E\^ = fi^ |k|^ /2m. Furthermore, 

1/2 

the distribution is sampled by a^^k = (^j,k + 1/2) /?j,k 
where Nj^y, = [exp((£;^-,k - l^j.o) /ksTo) - 1]""^ and 7/^-,k 
are the complex Gaussian random variables with mo- 
ments (r?j,k^i,k') = (^;,k^j*k') = ^ (^J,k^l,k') = 
(5kk'- Second, to simulate the thermal quench, the tem- 
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perature and chemical potential of the non-condensate 
band are altered to the new values T < Tq and /i > /ij^o- 
For convenience, we adopt the oscillator units in the nu- 
merical computations, where the len gth, tim e and energy 
are respectively scaled in units of ^/h/muJ^ uj~^ and huj. 

We first study the spinor BEC of ^^Rb, which has a 
Qs < 0. The total number of the modes are Ux^riy = 256 
and the energy cutoff is chosen at fixc^^yc = 128. The 
parameters are Q = 0.95, T = lOnK, /i = 8, and 
hjj/kBT = 0.03 for all ^j's. The time evolutions of 
the density profiles for ^j's are shown in Fig.l. During 
the evaporative cooling, the rotating condensates grow 
up and the emergent vortices start closely binding up 
and forming vortex-trimers in each . When the system 
reaches equilibrium, these vortex-trimers arrange them- 
selves into some interwoven lattice structures such that 
each vortex core of is filled up with particles of the 
rest two components. In other words, quantized vor- 
tices of either inter- or intraspecies, avoid to overlap 
with each other. The shapes of the trimer structures 
in all components are somewhat different. To charac- 
terize these structures, we calculate the incompressible 

kinetic energy per particle for each This can be inThlspkiorB^cIr^e'Afir « 4I5 x 10^^ 

done by writing = \'iij\exp{i(pj) and defining the 
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FIG. 3: (Color online) Snapshots of (a) |^-i|^ (b) |^o|^ 
(c) 1^1 1^; (d) 1^1^ for a quenched rotating spinor BEC of 
^^Na. When the system reaches equilibrium at T = 10 nK 
with Q = 0.8, /i = 25 (rightmost column), crystalline order of 
vortex-dimers is established in each ^ ^ . The particle numbers 
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the solenoidal and irrotational fields, where V • Z 
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0. The incompressible energy is defined by 

2 



(1/2) /^V 



which corresponds to the ki- 



netic energy of swirls in the superflows. Consequently, 
we find that sl^Nj = 19.48, 19.27, 20.05 for j = ±1,0, 
respectively. The spin textures are shown in Fig. 2 (a), 
where a hexagonal lattice is visualized. Furthermore, the 
spatial variation of the local spin moments parallel to 
the axis of rotation is plotted in Fig. 2(b). An enlarged 
perspective view of the 3-dimensional orientations of the 
local spins is shown in Fig. 2(c). Taking the island around 
the trap center for example, we see that the innermost 
spin points into the paper, while the others increasingly 
twist and bend in upwards direction. This vortex-like 
arrangement of magnetic moments is exactly the config- 
uration of a skyrmion [21-23 . Since the central spin in 
the skyrmion is perpendicular to the rotating plane, Sx 
and Sy must vanish therein, implying that the skyrmions 
must be centered at the regions of = 0, i.e., the cores 
of the vortex-trimer in ^o- The topological charge den- 
sity, a = s- {ds/dx X ds/dy) /47r [22 , is shown in Fig. 2 
(d), which exhibits a hexagonal lattice structure. Numer- 
ical integration over the primitive unit cell reveals that 
each skyrmion carries a topological charge Q = —1. We 
notice that the crystallization of skyrmions have recently 
been observed in various magnetic materials character- 
ized by the Dzyaloshinskii-Moriya interaction [22l [23] . 
Amazingly, the spin textures in these magnetic materials 
are highly resembling to those in Fig. 2 (a). 

For the case of gs > 0, we consider the spinor BECs 
of 2^Na. We set n = 0.8, T = 10 nK, /i = 25, and 
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FIG. 4: (Color online) (a) The equilibrium spin textures for 
the spinor BEC of ^^Na; (b) Distribution of Sz in the spm 
textures; (c) The orientations of the unit vecotr s (r) in the 
adjacent magnetic domains. The color of each arrow indicates 
the magnitude of 5'^; (d) The modulus of |^i - e-^^^^^^^^-i | , 
where dark shaded regions reveal the locations of HQ Vs. 



^Ij/ksT = 0.03 for all spin components. The total 
number of modes and energy cutoff remain the same as 
those for the case of ^^Rb. In Fig. 3, the time evolu- 
tions of the density profiles for are shown. The nu- 
cleation of the vortices in the current case are similar 
to that of the case of ^^Rb, except that a square lattice 
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of tightly bound vortex-dimers is formed in each com- 
ponent. Consequently, we find that 8^^-/Nj = 13.44, 
13.45, 12.54 for j = ±1,0, respectively. The equili- 
brated spin textures on the rotating plan are shown in 
Fig. 4(a), in which a mosaic of magnetic domains with 
staggered magnetization is created. In Fig. 4(b) and (c), 
we see that almost all spins in the magnetic domains lie 
in the x?/-plane. Spins belonging to the same domain 
align nearly unidirectionally. Note that the spins reverse 
their magnetization through a Bloch wall transition in a 
very narrow region near the boundary between two ad- 
jacent domains. These staggered magnetic domains act 
as the smoking gun of HQVs which have been predicted 
to exist in superfiuid ^He [2 and superconductors [24] . 
Considering the transformation, ^ G{0)R{n^x)^^ 
where G (0) = exp {iO) is a gauge transformation and 
R (n, x) = exp (^ix'^ • is a spin rotation through an 
angle x about the unit vector n, a HQV entails a spin 
rotation with x = '^r followed by a global phase change 
of TT in ^. Without loss of generality, we assume s (r) = 
cos (j) (r) Gx + sin (j) (r) e^. Numerically, we verify that the 
spin textures remain unchanged through the transfor- 
mation G (tt) ^ (s, tt) ^ = (e-2^^(^)^_i,^o,e2^^^''^^i)^. 

Upon requiring G (tt) ^ (s, tt) ^ =6*^"^^, it follows that 
g-2i0(r)^_^ = ^1 and e^^^W^^ ^ rj.^^ redun- 

dancy of the last equality implies that there are multiple 
solutions satisfying the criterion, |^i — e"^*'^^'^^^-! | = 0, 
or equivalently, |^-i| = |^i|. In Fig. 4(d), the modulus 
1^1 - e~^*'^'^'^^^_i| is plotted where the dark shaded ar- 
eas represent the core positions of HQVs, which form a 
square lattice apparently. Likewise, the HQVs are lo- 
calized at the cores of vortex-dimers in ^o- Our results 
are consistent with those obtained by means of dynami- 
cal creation [25^, where a lattice of HQVs can be created 
in a rotating optical trap when additional pulsed mag- 
netic trapping potentials are applied. Furthermore, the 
ground state of a rotating dipolar spinor BEC has been 
shown to have the same structure when the dipole-dipole 
interaction is small compared to the contact ones |26] . 



We note that an Q comparable to uu is needed to sta- 
bilize the crystalline orders of defects. Under such a fast 
rotation, the filling factor i.e., the ratio of the num- 
ber of atoms to the number of vortices, can have a value 
of few hundreds for each component, as shown in Fig.l 
and Fig. 3. According to the criterion in [27 , the system 
enters the mean-field quantum Hall regime, in which the 
mean-field theory still applies yet the state of the system 
can be well described in the lowest Landau level approxi- 
mation. When Q is not sufficiently large, the crystalliza- 
tion does not arise albeit some few topological defects 
may be readily created in the condensate. For example, 
in our simulations with Q = 0.3, we find only a few HQVs 
nucleating in the spinor BEC of ^^Na during the rapid ro- 
tational evaporative cooling. The situation is somewhat 
different in the case of ^^Rb, where Mermin-Ho vortices, 
rather than the skyrmions, are created in the condensate. 
Furthermore, to see whether the crystallization is robust 
against the thermal ffuctuations arising from the growth 
of condensate, we have assumed a range of final equilib- 
rium temperatures in the simulations. We find that both 
crystalline orders remain intact for temperatures lower 
than 50nK. At higher temperatures, however, the crys- 
tallization is thwarted by the fiuctuations of spin textures 
so that the lattice becomes disordered and starts melting 
when the system approaches the critical regime. 

In summary, we have investigated the non-equilibrium 
dynamics of spin-1 BECs during the rapid rotational 
evaporative cooling. Crystallization of skyrmions and 
HQVs is predicted to arise in the spinor BEC of ^^Rb 
and ^^Na, respectively. To resolve the spatial magneti- 
zation of the crystallized topological defects, the images 
have to be taken in situ, and this can be achieved basi- 
cally by using the polarization-dependent phase-contrast 
technique [28 . 
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